//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   testPsi_Cell_Field.cc
 * \author Jeremy Roberts
 * \date   Aug 17, 2011
 * \brief  Test of Psi_Cell_Field
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev:: 130                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-09-07 02:55:28 +0000 (Wed, 07 Sep 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>
#include <cmath>

#include "test/Testing.hh"

#include "tran/Psi_Cell_Field.hh"
#include "tran/CartMesh.hh"
#include "tran/Traits.hh"
#include "tran/Definitions.hh"
#include "../utilities/Constants.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"
#include "../utilities/SP.hh"
// Equations
#include "tran/DD_1D_Equations.hh"
#include "tran/DD_2D_Equations.hh"
#include "tran/DD_3D_Equations.hh"

using namespace std;
using namespace slabtran;
using namespace util;

int test1d_1();
int test2d_1();
int test3d_1();

int main(int argc, char *argv[])
{
  Require(argc==2);
  int test = atoi(argv[1]);
  cout << " ------ testPsi_Cell_Field ------ " << endl;
  int test_value = 0;
  int number_of_tests = 4;
  if (test == 1 || test == number_of_tests) test_value += test1d_1();
  if (test == 2 || test == number_of_tests) test_value += test2d_1();
  if (test == 3 || test == number_of_tests) test_value += test3d_1();

  if (test > number_of_tests || test < 1) cout << " test " << test
      << " not available. " << endl;
  if (test_value) cout << " failed with value = " << test_value << endl;
  else
    cout << " passed." << endl;
  return test_value;
  return 0;
}

int test1d_1()
{

  // make a simple mesh
  double width = 10.0;
  int count = 10;
  Mesh<_1D>::SP_mesh mesh;
  mesh = new Mesh<_1D>(width, count);

  // initialize the numbers of angles and octants
  int number_angles_oct = 2;
  int number_octants = 2;
  int number_angles = number_angles_oct * number_octants;

  Psi_Cell_Field<DD_1D_Equations,_1D>::SP_psi_cell_field flux;
  flux = new Psi_Cell_Field<DD_1D_Equations,_1D>(mesh, number_angles, 0.0);

  //cout << " flux size = " << flux->size() << endl;
  int c = 0; // cardinal index counter
  for (int o = 0; o < number_octants; o++)
    for (int n = 0; n < number_angles_oct; n++)
       for (int i = 0; i < count; i++, c++)
         (*flux)[c] += 100.0 * i + 10.0 * o + 1.0 * n * n;

  TEST(c == flux->size());
  TEST(number_angles == flux->number_angles());

  // verify
  for (int i = 0; i < count; i++)
    for (int o = 0; o < number_octants; o++)
      for (int n = 0; n < number_angles_oct; n++)
      {
        int angle = o * number_angles_oct + n;
        TEST((*flux)(i, 0, angle) == 100.0 * i + 10.0 * o + 1.0 * n * n);
      }

  cout << " ... passed 1d psi cell field test. " << endl;
  return 0;
}

int test2d_1()
{

  // make a simple mesh
  double width = 10.0;
  int count = 10;
  Mesh<_2D>::SP_mesh mesh;
  mesh = new Mesh<_2D>(width, count);

  // initialize number of angles
  int number_angles_oct = 4;
  int number_octants = 4;
  int number_angles = number_angles_oct * number_octants;

  // total size should be (10*10)*7 = 700
  Psi_Cell_Field<DD_2D_Equations,_2D> flux(mesh, number_angles, 0.0);

  /*!
   * The following is the suggested way to update the moments
   * within the equation solver (e.g. \ref DD_2D_Equations).
   * The programmer implementing a 2d solver is \em obligated
   * to update only those moments valid for 2d, which means
   * \em not trying to update \f$\phi^l_0\f$ for \f$l>0\f$. This
   * condition is easily estatblished if the 2d moments interface
   * is used, as it requires the caller uses valid indices; hence,
   * the programmer will be quickly alerted to any errors.  So,
   * use the 2d interface!
   *
   */
  int c = 0;
  for (int o = 0; o < number_octants; o++)
  {
    for (int n = 0; n < number_angles_oct; n++)
    {
      for (int j = 0; j < count; j++)
      {
        for (int i = 0; i < count; i++)
        {
          flux[c++] += 10.0 * i + 100.0 * j + 1000 * o + 10000.0 * n;
        }
      }
    }
  }
  TEST(c == flux.size());
  TEST(number_angles == flux.number_angles());

  // verify
  for (int j = 0; j < count; j++)
  {
    for (int i = 0; i < count; i++)
    {
      for (int o = 0; o < number_octants; o++)
      {
        for (int n = 0; n < number_angles_oct; n++)
        {
          int angle = o * number_angles_oct + n;
          TEST(flux(i, j, 0, angle) ==
               10.0 * i + 100.0 * j + 1000 * o + 10000.0 * n);
        }
      }
    }
  }
  cout << " ... passed 2d psi cell field test. " << endl;

  return 0;
}

int test3d_1()
{

  // make a simple mesh
  double width = 10.0;
  int count = 10;
  Base_Mesh::SP_mesh mesh; // use base pointer
  mesh = new Mesh<_3D>(width, count);

  // initialize numbers of angles and octants
  int number_angles_oct = 4;
  int number_octants = 8;
  int number_angles = number_angles_oct * number_octants;

  Psi_Cell_Field<DD_3D_Equations,_3D> flux(mesh, number_angles, 0.0);

  int c = 0;

  for (int o = 0; o < number_octants; o++)
  {
    for (int n = 0; n < number_angles_oct; n++)
    {
      for (int k = 0; k < count; k++)
      {
        for (int j = 0; j < count; j++)
        {
          for (int i = 0; i < count; i++)
          {
            flux[c++] += 1. * i + 10. * j + 100. * k + 1000. * o + 10000. * n;
          }
        }
      }
    }
  }

  TEST(c == flux.size());
  TEST(number_angles == flux.number_angles());

  // verify
  for (int k = 0; k < count; k++)
  {
    for (int j = 0; j < count; j++)
    {
      for (int i = 0; i < count; i++)
      {
        for (int o = 0; o < number_octants; o++)
        {
          for (int n = 0; n < number_angles_oct; n++)
          {
            int angle = o * number_angles_oct + n;
            TEST(flux(i, j, k, 0, angle)
                 == 1.0 * i + 10.0 * j + 100.0 * k + 1000 * o + 10000.0 * n);
          }
        }
      }
    }
  }
  cout << " ... passed 3d psi cell field test. " << endl;

  return 0;
}

//---------------------------------------------------------------------------//
//              end of testPsi_Cell_Field.cc
//---------------------------------------------------------------------------//
